Who Am I?

Hi! My name is Dr. Elizabeth Sweeney. I am currently an Assistant Professor in the Biostatistics division at Weill Cornell and today marks my 4th day on the job! The journey that led me to where I am now was a little winding…

My research interests are in working with electronic health records (EHR) data and structural magnetic resonance imaging (sMRI) data. When I’m not doing research I like hiking and riding my bike.

A few fun facts about me:

Hiking in the Catskills in the freezing cold to get my winter peaks.

Hiking in the Catskills in the freezing cold to get my winter peaks.

What is Neuroconductor?

Here is a link to the Neuroconductor website: https://neuroconductor.org

Neuroconductor is an open-source platform for rapid testing and dissemination of reproducible computational imaging software. The goals of the project are to:

Based on the programming language R, Neuroconductor started with 51 inter-operable packages that cover multiple areas of imaging including visualization, data processing and storage, and statistical inference. Neuroconductor accepts new R package submissions, which are subject to a formal review and continuous automated testing.

In short, all of the good packages for doing brain image analysis in R are hosted on Neuroconductor!

The Core Team of Neuroconductor?

Neuroconductor was developed at Johns Hopkins. The Core Team of Neuroconductor includes:

John Muschelli, R package guru Adi Gherman, Lead developer Brian Caffo, Professor Ciprian Crainiceanu, Professor

How to Install Neuroconductor

You can install neuroconductor using the commands below. If you need more assistance please visit the instructions at the link below:

https://neuroconductor.org/tutorials/install

##installation of neruoconductor requires the installation of the package 'devtools' 
install.packages('devtools', repos = "http://cran.us.r-project.org")
## 
##   There is a binary version available but the source version is
##   later:
##          binary source needs_compilation
## devtools 1.13.6  2.0.2             FALSE
## installing the source package 'devtools'
##code to install nueroconductor 
source("https://neuroconductor.org/neurocLite.R")
## Bioconductor version 3.6 (BiocInstaller 1.28.0), ?biocLite for help
## A new version of Bioconductor is available after installing the most
##   recent version of R; see http://bioconductor.org/install
## BioC_mirror: https://bioconductor.org
## Using Bioconductor 3.6 (BiocInstaller 1.28.0), R 3.4.4 (2018-03-15).
## Using neurocLite version: # neurocInstall package version: 0.11.2

How to Install a Package with Neuroconductor

Let’s use neuroconductor to install all of the packages we will need for this tutorial. Note that I have commented out the ‘neuro_install’ command as this did not play well with an Rmarkdown document :( When you are working on your own please remove the commenting out and install all of these packages from Neuroconductor.

#neuro_install('oro.nifti')
library(oro.nifti)
## oro.nifti 0.9.11
## 
## Attaching package: 'oro.nifti'
## The following object is masked from 'package:dplyr':
## 
##     slice
#neuro_install('fslr')
library(fslr) 
## Loading required package: neurobase
#neuro_install('WhiteStripe')
library(WhiteStripe)
#neuro_install('oasis')
library(oasis)

Neuroimaging Informatics Technology Initiative (NIfTI) File Format

Much of neuroimaging data is stored in a NIfTI file format. These files have a .nii or .nii.gz extension.

The Package oro.nifti Works with NIfTI Files

The R package oro.nifti:

We will use a subject from the Kirby 21, an open source multi-modal sMRI reproducibility study with 21 healthy subjects (www.nitrc.org/projects/multimodal) to show what we can do with the oro.nifti package:

We start by using the ‘readNIfTI’ function to read in the baseline T1-weighted sMRI from the Kirby 21 study. A T1-weighted sMRI is one modality of sMRI; we will see more modalities later in the tutorial.

MPRAGE_base <- readNIfTI('SUBJ0001-01-MPRAGE.nii.gz', reorient=FALSE)

The image is stored as an 3-dimensional array. Each element in the array is referred to as a ‘voxel’. Each voxel takes on a numeric value.

dim(MPRAGE_base)
## [1] 170 256 256

We can get header information from the image. This tells us information about how the image was acquired on the MRI scanner and how it is stored.

MPRAGE_base
## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 16 (FLOAT32)
##   Bits per Pixel  : 32
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 1 (Scanner_Anat)
##   Sform Code      : 0 (Unknown)
##   Dimension       : 170 x 256 x 256
##   Pixel Dimension : 1.2 x 1 x 1
##   Voxel Units     : mm
##   Time Units      : sec

A nifti R object is an S4 object. Let’s check out the slot names for the S4 object.

slotNames(MPRAGE_base)
##  [1] ".Data"          "sizeof_hdr"     "data_type"      "db_name"       
##  [5] "extents"        "session_error"  "regular"        "dim_info"      
##  [9] "dim_"           "intent_p1"      "intent_p2"      "intent_p3"     
## [13] "intent_code"    "datatype"       "bitpix"         "slice_start"   
## [17] "pixdim"         "vox_offset"     "scl_slope"      "scl_inter"     
## [21] "slice_end"      "slice_code"     "xyzt_units"     "cal_max"       
## [25] "cal_min"        "slice_duration" "toffset"        "glmax"         
## [29] "glmin"          "descrip"        "aux_file"       "qform_code"    
## [33] "sform_code"     "quatern_b"      "quatern_c"      "quatern_d"     
## [37] "qoffset_x"      "qoffset_y"      "qoffset_z"      "srow_x"        
## [41] "srow_y"         "srow_z"         "intent_name"    "magic"         
## [45] "extender"       "reoriented"

Visualizing images with the base R image function

We can visualize a single slice of the image in different orientations using the ‘image’ function in R.

##axial slice##
image(MPRAGE_base[,,128])

##coronal slice##
image(MPRAGE_base[,128,], col = rainbow(12))

##sagittal slice##
image(MPRAGE_base[85,,], col = topo.colors(12))

Visualizing images using the orthographic function (from the oro.nifti package)

These graphics that we made with image are just okay… Let’s use the ‘orthographic’ function from the oro.nifti package to make better looking graphics.

orthographic(MPRAGE_base)

The default of the orthographic function is to take the middle slice in each direction (coronal, sagital and axial direction). We can also pass in the slices which we would like to use:

orthographic(MPRAGE_base, xyz = c(90, 100, 15))

Make a histogram of the data

The sMRI image is just a 3-dimensional array of numbers. Let’s look at a summary of the data in the image.

summary(MPRAGE_base)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##       0       0    3416  143789  256204 2408322

We can also use regular plot functions to do things like make histograms of the data.

ggplot.MPRAGE.object <- data.frame(intensities = c(MPRAGE_base)) 

ggplot(ggplot.MPRAGE.object, aes(x = intensities)) + 
  geom_histogram() 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Most of the intensities in the background of the image are set to 0, we can remove these to get a better histogram.

ggplot.MPRAGE.object.no.zeros <- ggplot.MPRAGE.object  %>% 
  filter(intensities > 0)

dim(ggplot.MPRAGE.object)
## [1] 11141120        1
dim(ggplot.MPRAGE.object.no.zeros)
## [1] 5981646       1
ggplot(ggplot.MPRAGE.object.no.zeros, aes(x = intensities)) + 
  geom_histogram() 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Motivation for image pre-processing

Let’s read in a follow-up image from the Kirby 21 dataset for this same subject.

MPRAGE_follow <- readNIfTI('SUBJ0001-02-MPRAGE.nii.gz', reorient = FALSE)

We can view the baseline and follow-up images side by side using the ‘double_ortho’ function from fslr

double_ortho(MPRAGE_base, MPRAGE_follow)

We might be curious what the difference between these two images is. We can subtract the two images and see:

MPRAGE_diff <- MPRAGE_base - MPRAGE_follow

Note that the subtraction of the two images is still a nifti R object

MPRAGE_diff 
## NIfTI-1 format
##   Type            : nifti
##   Data Type       : 16 (FLOAT32)
##   Bits per Pixel  : 32
##   Slice Code      : 0 (Unknown)
##   Intent Code     : 0 (None)
##   Qform Code      : 1 (Scanner_Anat)
##   Sform Code      : 0 (Unknown)
##   Dimension       : 170 x 256 x 256
##   Pixel Dimension : 1.2 x 1 x 1
##   Voxel Units     : mm
##   Time Units      : sec

Visualize the difference of the two images:

orthographic(MPRAGE_diff )

The images are not exactly aligned – which creates this ghosting effect. To assess true differences we need to do some image pre-processing.

The fslr package for preprocessing data

Here is a link to the information about the fslr package: https://neuroconductor.org/package/fslr

FSL is useful, open-source, scriptable software for neuroimaging analysis. One of the main problems with using FSL is that it requires coding in bash. The solution is to use the fslr package. The fslr package ports many of the main functions of FSL into R. A Big Disclaimer: May not work on Windows operating systems!

To use fslr you first need to install FSL. Here are the steps to install FSL and start using it with fslr (note that we have already installed the fslr package in this tutorial).

The steps of image processing are outlined below:

Inhomogenity correction

The first step that we will do is to remove any image Inhomogeneity. Image inhomogeneity is background image intensity variation that is caused by imperfections in imaging devices. This is the first step of our image pre-processing pipeline

To use FSL, I first need to set my options to my FSL path

options(fsl.path= "/usr/local/fsl")

Now I am ready to use the ‘fsl_biascorrect’ function to do the inhomogeneity correction. Note that we won’t do this during the tutorial as it takes some time to run the algorithm!

## MPRAGE_base_bias_corrected <- fsl_biascorrect(MPRAGE_base) 

## writeNIfTI(MPRAGE_base_bias_corrected, filename = 
##            'SUBJ0001-01-MPRAGE_bias_corr', verbose = TRUE, gzipped = TRUE)

## MPRAGE_follow_bias_corrected <- fsl_biascorrect(MPRAGE_follow) 

## writeNIfTI(MPRAGE_follow_bias_corrected, filename = 
##            'SUBJ0001-02-MPRAGE_bias_corr', verbose = TRUE, gzipped = TRUE)

We will now read in the inhomogeneity (or bias corrected) image. We will subtract it from the original image to see the removed inhomogeneity.

## read in the bias corrected baseline image 
MPRAGE_base_bias_corrected <- readNIfTI('SUBJ0001-01-MPRAGE_bias_corr.nii.gz', reorient=FALSE)

##take a difference between the bias corrected image and the original image 
bias_diff <- MPRAGE_base_bias_corrected  - MPRAGE_base

##visualize the difference 
orthographic(bias_diff)

Skull stripping

Skull stripping removes unwanted tissue (the skull, neck, etc.) and leaves behind only the brain.

We now perform skull stripping with fslbet

MPRAGE_base_bias_corrected_stripped <- fslbet(MPRAGE_base_bias_corrected) 
## Warning in get.fsloutput(): Can't find FSLOUTPUTTYPE, setting to NIFTI_GZ
## FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/bet2 "/private/var/folders/hx/5nm8nvg93pq252g_q5z50r_m0000gn/T/RtmpncSWyv/file12ee83edd72cb.nii.gz" "/var/folders/hx/5nm8nvg93pq252g_q5z50r_m0000gn/T//RtmpncSWyv/file12ee85ac8f34d"

We can see some of the inner workings of the fslr package by looking more closely at the fslbet command. We can see how fslr uses system calls to call FSL.

fslbet
## function (infile, outfile = NULL, retimg = TRUE, reorient = FALSE, 
##     intern = FALSE, opts = "", betcmd = c("bet2", "bet"), verbose = TRUE, 
##     ...) 
## {
##     betcmd = match.arg(betcmd)
##     cmd <- get.fsl()
##     outfile = check_outfile(outfile = outfile, retimg = retimg, 
##         fileext = "")
##     infile = checkimg(infile, ...)
##     outfile = checkimg(outfile, ...)
##     outfile = nii.stub(outfile, ...)
##     cmd <- paste0(cmd, sprintf("%s \"%s\" \"%s\" %s", betcmd, 
##         infile, outfile, opts))
##     if (verbose) {
##         message(cmd, "\n")
##     }
##     res = system(cmd, intern = intern)
##     ext = get.imgext()
##     outfile = paste0(outfile, ext)
##     if (retimg) {
##         img = readnii(outfile, reorient = reorient, ...)
##         return(img)
##     }
##     return(res)
## }
## <environment: namespace:fslr>

Let’s return to our skull stripped image. We can plot the stripped image next to the original image using the ‘double_ortho’ function:

double_ortho(MPRAGE_base, MPRAGE_base_bias_corrected_stripped)

We can also make a ‘brain mask’ from the skull stripped image

bet_mask <- niftiarr(MPRAGE_base_bias_corrected_stripped, 1)
is_in_mask = MPRAGE_base_bias_corrected_stripped>0
bet_mask[!is_in_mask]<-0

We can plot the brain mask.

orthographic(bet_mask)

We can now make a histogram of just the data after skull stripping.

##add the brain mask to our dataframe 'ggplot.MPRAGE.object'
ggplot.MPRAGE.object$bet_mask <- c(bet_mask)

##filter the data frame for voxels that take a value 1 in the brain mask 
ggplot.MPRAGE.object.masked <- ggplot.MPRAGE.object %>% 
  filter(bet_mask == 1)

## plot the histogram of the data 
ggplot(ggplot.MPRAGE.object.masked, aes(x = intensities)) + 
  geom_histogram() 
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Or we can make a density plot, where we can see two distinct peaks one for gray matter and one for white matter!

ggplot(ggplot.MPRAGE.object.masked, aes(x = intensities)) + 
  geom_density()

Now we can write the brain mask using the ‘writeNIfTI’ function from the oro.nifti package:

writeNIfTI(bet_mask, filename = 
             "brain_mask", verbose = TRUE, gzipped = TRUE)
##   writing data at byte = 352

Image Registration

We now register the two images. There are many different registration algorithms that you can use with fsl (both linear and non-linear). We will use a linear registration with flirt. This registration has 6 degrees of freedom (roll, pitch and yaw). This is appropriate because we have a baseline and follow-up image for a patient and therefore do not need to any non-linear transformations.

MPRAGE_follow_bias_corrected <- readNIfTI('SUBJ0001-02-MPRAGE_bias_corr.nii.gz', reorient=FALSE)


MPRAGE_follow_reg_base <- flirt(MPRAGE_follow_bias_corrected, 
                                MPRAGE_base_bias_corrected, 
                                dof = 6, 
                                retimg = TRUE, 
                                reorient = FALSE)
## FSLDIR='/usr/local/fsl'; PATH=${FSLDIR}/bin:${PATH};export PATH FSLDIR; sh "${FSLDIR}/etc/fslconf/fsl.sh"; FSLOUTPUTTYPE=NIFTI_GZ; export FSLOUTPUTTYPE; ${FSLDIR}/bin/flirt -in "/private/var/folders/hx/5nm8nvg93pq252g_q5z50r_m0000gn/T/RtmpncSWyv/file12ee83521a0.nii.gz" -ref "/private/var/folders/hx/5nm8nvg93pq252g_q5z50r_m0000gn/T/RtmpncSWyv/file12ee820328f7b.nii.gz" -out "/var/folders/hx/5nm8nvg93pq252g_q5z50r_m0000gn/T//RtmpncSWyv/file12ee86f6d88df" -dof 6 -omat "/var/folders/hx/5nm8nvg93pq252g_q5z50r_m0000gn/T//RtmpncSWyv/file12ee85769e4c0.mat"

Now we can take the difference of the baseline and follow-up images:

MPRAGE_reg_diff <- MPRAGE_base_bias_corrected - MPRAGE_follow_reg_base

And visualize this versus our original differencing without registration:

double_ortho(MPRAGE_diff, MPRAGE_reg_diff)

Intensity Normalization

The next step is to intensity normalize the images. sMRI images do not have units, so one method of normalizing these images is to put them into z-score units of the normal appearing white matter. We do this with the R package WhiteStripe. A really nice tutorial for WhiteStripe can be found here: https://cran.r-project.org/web/packages/WhiteStripe/vignettes/Running_WhiteStripe.html/

##
ws = whitestripe(img = MPRAGE_base_bias_corrected_stripped, 
                 type = "T1", 
                 stripped = TRUE)
## Making Image VOI
## Making T1 Histogram
## Getting T1 Modes
## Smoothing Histogram
## Smoothing Derivative
## Quantile T1 VOI
norm = whitestripe_norm(img =  MPRAGE_base_bias_corrected_stripped, 
                        indices = ws$whitestripe.ind)

dim(norm) 
## [1] 170 256 256

When we look at the image we don’t see much difference from the original image (because this is a z-score normalization, the image is just shifted and scaled).

orthographic(norm)

Where we will see a difference is in the densities of the intensities in the image.

ggplot.MPRAGE.object$normalized <- c(norm)


ggplot.MPRAGE.object.norm.masked <- ggplot.MPRAGE.object %>% 
  filter(bet_mask == 1)

Here are the original intensities

ggplot(ggplot.MPRAGE.object.norm.masked, aes(x = intensities)) + 
  geom_density()

And here are the intensities after white stripe normalization

ggplot(ggplot.MPRAGE.object.norm.masked, aes(x = normalized)) + 
  geom_density()

Resources!

Want to learn more about image pre-processing and visualization?

Multiple Sclerosis Lesion Segmentation

sMRI Imaging Modalities

So far, we have been working with T1-weighted images from the Kirby 21 study. There are different modalities for sMRI images that produce different images. When a person goes to have an sMRI study taken they will often have multiple sMRIs with different modalities taken of their brain.

We will now switch to a different dataset from a patient with multiple sclerosis. The data can be found here https://github.com/muschellij2/imaging_in_r/tree/gh-pages/output/orig

Let’s look at a Fluid Attenuated Inversion Recovery (FLAIR) image from a patient with Multiple Sclerosis (MS):

FLAIR <- readNIfTI('training01_01_flair.nii.gz', reorient=FALSE)

The hyperintense areas of the image are brain lesions.

orthographic(FLAIR)

We can look at different slices to see more brain lesions.

orthographic(FLAIR, xyz = c(100, 100, 80))

Now let’s look at the T2-weighted image

T2 <- readNIfTI('training01_01_t2.nii.gz', reorient=FALSE)

Again, lesions are hyperintense on the T2-weighted image

orthographic(T2)

Now let’s look at the Proton Density (PD) image

PD <- readNIfTI('training01_01_pd.nii.gz', reorient=FALSE)

Lesions are also hyperintense on the T2-weighted image.

orthographic(PD)

And we can also look at the T1-weighted image

T1 <- readNIfTI('training01_01_mprage_n4.nii.gz', reorient=FALSE)

Lesions are hypointense on the T1-weighted image.

orthographic(T1)

Next we use the ‘oasis_predict’ function to make a map of where the lesions in the brain are. This is an algorithm that we developed which takes a logistic regression model and trains on gold standard lesion Segmentations made by a radiologist. Different features for the logistic regression model are generated using the FLAIR, PD, T2, and T1 images.

Link for the OASIS paper: https://www.sciencedirect.com/science/article/pii/S2213158213000235

The ‘oasis_predict’ function takes some time to run so again we comment this out and I have provided the output from the ‘oasis_predict’ function.

##oasis_map <- oasis_predict(flair = FLAIR ,
##                           t1 = T1,
##                           t2 = T2,
##                           pd = PD)


## writeNIfTI(oasis_map$oasis_map, 
## filename = 'oasis', verbose = TRUE, gzipped = TRUE)

We will read in the output from the oasis_predict function. We now plot the probability map of lesion presence that is generated by oasis_predict.

oasis_map <- readNIfTI('oasis.nii.gz', reorient=FALSE)

orthographic(oasis_map)

We can binarise the probability map to create a binary map of lesion presence.

oasis_map_binary <- oasis_map
oasis_map_binary[oasis_map_binary >= .5] <- 1
oasis_map_binary[oasis_map_binary < 1] <- NA

And we can plot this over the T1-weighted image to see if we located the lesions.

orthographic(T1, oasis_map_binary, col.y = "orange")

Let’s check another set of slices to assess performance further!

orthographic(FLAIR, xyz = c(100, 100, 80))

orthographic(T1, oasis_map_binary, col.y = "orange", xyz = c(100, 100, 80))